University of Iowa

College of Liberal Arts and Sciences

Basic Data Wrangling and Data Visualization in R

Introduction

In this class I will

• Briefly outline the history of R.
• Using some examples briefly show how to do data wrangling and visualize data in R.

Materials for this class are available on GitHub at https://github.com/ltierney/SIBS-WV-2021.git.

• You can access it as an RStudio project by following the menu selection File > New Project > Version Control > Git and specifying this URL.

• You can use the git command line client with

  git clone https://github.com/ltierney/SIBS-WV-2021.git

Materials for our Data Visualization and Data Technologies course are available at http://www.stat.uiowa.edu/~luke/classes/STAT4580-2021/

Some tools I will be using:

• The RStudio IDE.
• Many features from the basic R distribution.
• Some tools from the tidyverse.
• The ggplot package based on the Grammar of Graphics framework.

Most of the packages are loaded by loading the tidyverse package:

library(tidyverse)

Useful references:

Hadley Wickham and Garrett Grolemund (2016), R for Data Science, O'Reilly.

Claus O. Wilke (2019), Fundamentals of Data Visualization, O'Reilly.

Kieran Healy (2018) Data Visualization: A practical introduction, Princeton

Rafael A. Irizarry (2019), Introduction to Data Science: Data Analysis and Prediction Algorithms with R, Chapman & Hall/CRC. (Book source on GitHub)

Ask questions any time!

The R Language

R is a language for data analysis and graphics.

• R was originally developed by Robert Gentleman and Ross Ihaka in the early 1990's for a Macintosh computer lab at U. of Auckland, New Zealand.

• R is based on the S language developed by John Chambers and others at Bell Labs.

R is an Open Source project.

• Since 1997 R is developed and maintained by the R-core group, with around 20 members located in maor than 10 different countries.

• R is widely used in the field of statistics and beyond, especially in university environments.

• R has become the primary framework for developing and making available new statistical methodology.

• Many (now over 17,000) extension packages are available through CRAN or similar repositories.

Working with R

R is designed for interactive data exploration.

• Interaction is through a read-eval-print loop (REPL).
• This is also called a command line interface (CLI).

All computations are specified in the R language.

• Even for simple tasks you need to know a little of the language.
• After learning to do simple tasks you know some of the language.

The language is used to

• prepare data for analysis;
• specify individual analyses;
• program repeated or similar analyses;
• program new methods of analysis.

Specifying these tasks in a language supports reproducible research.

The R language operates on vectors and arrays.

Commonly used data types are:

• integer and numeric vectors;
• logical vectors;
• character vectors;
• factors.

All basic vector types support missing (NA) values.

Arithmetic operations are vectorized to operate element-wise on vectors.

Data vectors are usually combined into table-like objects called data frames.

The Data Analysis Process

A figure that shows the steps usually involved in a data analysis project:

These steps are often repeated many times, so it is important to make your work reproducible.

Reproducible Data Analysis

Making your work reproducible:

• Save you work in a text file or notebook.

• Track changes to your files with a version control system like git.

• Use a system like Rmarkdown to prepare your reports.

This allows you to re-create your report when data changes (as it often will!).

A good resource for setting up your tools to support this is Happy Git and GitHub for the useR.

Some Examples

Working with research data a first step is usually to read and clean the data.

We'll put that off for a little while and work with some data sets made available in R packages.

Data sets available in R packages include:

• many classic data sets;
• newer, often larger, data sets useful for learning;
• current data obtained by querying web APIs.

Old Faithful Eruptions

A simple classic data set is the geyser data frame available in package MASS

data(geyser, package = "MASS")
dim(geyser)
## [1] 299   2
head(geyser)
##   waiting duration
## 1      80 4.016667
## 2      71 2.150000
## 3      57 4.000000
## 4      80 4.000000
## 5      75 4.000000
## 6      77 2.000000

head and tail return the first and last few rows of a data frame.

They are useful for quick sanity checks.

The rows represent measurements recorded for eruptions of the Old Faithful geyser in Yellowstone National Park, Wyoming.

The variables are:

• waiting: the time in minutes since the previous eruption;
• duration: the duration of the eruption.

The durations have a bimodal distribution:

ggplot(geyser) +
    geom_histogram(aes(x = duration), bins = 15, color = "black", fill = "grey")

A basic template for creating a plot with ggplot:

ggplot(data = <DATA>) + <GEOM>(mapping = aes(<MAPPINGS>))

An interesting question is whether the duration can be used to predict when the next eruption will occur.

A plot of the previous duration against the waiting time to the current eruption:

ggplot(geyser) + geom_point(aes(x = lag(duration), y = waiting))
## Warning: Removed 1 rows containing missing values (geom_point).

It looks like a useful rule would be to expect a shorter waiting time after a shorter eruption.

An interesting feature: Many durations are recorded as 2 or 4 minutes. This can also be seen in a histogram with many small bins:

p <- ggplot(geyser) +
    geom_histogram(aes(x = duration, y = stat(density)),
                   fill = "grey", color = "black", bins = 50)
p

ggplot produces a plot object. Drawing only happens when the object is printed.

Does this rounding matter?

• For many analyses it probably doesn't.
• It might if you wanted to fit normal distributions to the two groups.

Taking 3 minutes as the divide between short and long durations we can compute the means and standard deviations as

d <- geyser$duration
d_short <- d[d < 3]
d_long <- d[d >= 3]
mean(d_short)
## [1] 1.980317
sd(d_short)
## [1] 0.2779829
mean(d_long)
## [1] 4.262113
sd(d_long)
## [1] 0.3937525
mean(d >= 3)
## [1] 0.6488294

An approach that scales better is to compute group summaries using tools from the dplyr tidyverse package.

First, add a type variable:

geyser <- mutate(geyser, type = ifelse(duration < 3, "short", "long"))

The summaries can then be computed as

sgd <- summarize(group_by(geyser, type),
                 mean = mean(duration),
                 sd = sd(duration),
                 n = n())
(sgd <- mutate(sgd, prop = n / sum(n)))
## # A tibble: 2 x 5
##   type   mean    sd     n  prop
##   <chr> <dbl> <dbl> <int> <dbl>
## 1 long   4.26 0.394   194 0.649
## 2 short  1.98 0.278   105 0.351

This computation can also be written using the forward pipe operator %>%:

sgd <-
    group_by(geyser, type) %>%
    summarize(mean = mean(duration),
              sd = sd(duration),
              n = n()) %>%
    ungroup() %>%
    mutate(prop = n / sum(n))

The pipe operator allows a sequence of operations to be chained together.

The left-hand operation is passed implicitly as the first argument to the function called on the right.

One way to show the superimposed normal densities:

f1 <- function(x)
    sgd$prop[1] * dnorm(x, sgd$mean[1], sgd$sd[1])
f2 <- function(x)
    sgd$prop[2] * dnorm(x, sgd$mean[2], sgd$sd[2])
p <- p +
    stat_function(color = "red", fun = f1) +
    stat_function(color = "blue", fun = f2)
p

A ggplot can consist of several layers.

The means and standard deviations are affected by the rounding. Summaries that omit values equal to 2 or 4 minutes can be computed as

geyser2 <- filter(geyser, duration != 2, duration != 4)
sgd2 <-
    group_by(geyser2, type) %>%
    summarize(mean = mean(duration),
              sd = sd(duration),
              n = n()) %>%
    ungroup() %>%
    mutate(prop = n / sum(n))
sgd2
## # A tibble: 2 x 5
##   type   mean    sd     n  prop
##   <chr> <dbl> <dbl> <int> <dbl>
## 1 long   4.36 0.422   141 0.632
## 2 short  1.97 0.315    82 0.368

summarize, group_by, and mutate are from the dplyr package that implements a grammar of data manipulation.

A plot showing curves computed both ways:

f1_2 <- function(x)
    sgd2$prop[1] * dnorm(x, sgd2$mean[1], sgd2$sd[1])
f2_2 <- function(x)
    sgd2$prop[2] * dnorm(x, sgd2$mean[2], sgd2$sd[2])
p <- p +
    stat_function(color = "red",
                  linetype = 2,
                  fun = f1_2) +
    stat_function(color = "blue",
                  linetype = 2,
                  fun = f2_2)
p

Minnesota Barley Yields

A classic data set: Total yield in bushels per acre for 10 varieties at 6 sites in Minnesota in each of two years, 1931 and 1932.

The raw data:

data(barley, package = "lattice")
head(barley)
##      yield   variety year            site
## 1 27.00000 Manchuria 1931 University Farm
## 2 48.86667 Manchuria 1931          Waseca
## 3 27.43334 Manchuria 1931          Morris
## 4 39.93333 Manchuria 1931       Crookston
## 5 32.96667 Manchuria 1931    Grand Rapids
## 6 28.96667 Manchuria 1931          Duluth

Some initial plots:

p1 <- ggplot(barley) + geom_point(aes(x = yield, y = variety))
p2 <- ggplot(barley) + geom_point(aes(x = yield, y = site))
cowplot::plot_grid(p1, p2)

Using color to separate yields in the two years:

p1 <- ggplot(barley) + geom_point(aes(x = yield, y = variety, color = year))
p2 <- ggplot(barley) + geom_point(aes(x = yield, y = site, color = year))
cowplot::plot_grid(p1, p2)

Can we also show site using symbol shape?

ggplot(barley) +
    geom_point(aes(x = yield, y = variety, color = year, shape = site))

There is a lot of interference between shape and color.

Possible improvements:

• jittering;
• larger points.

ggplot(barley) +
    geom_point(aes(x = yield, y = variety, color = year, shape = site),
               position = position_jitter(height = 0.15, width = 0),
               size = 2)

Another approach: faceting to produce small multiples.

ggplot(barley) +
    geom_point(aes(x = yield, y = variety, color = year)) +
    facet_wrap(~site)

Focusing on summaries can help. Bar charts are sometimes used for summaries, but dot plots are usually a better choice.

barley_site_year <-
    group_by(barley, site, year) %>%
    summarize(yield = mean(yield)) %>%
    ungroup()
## `summarise()` has grouped output by 'site'. You can override using the `.groups` argument.
p1 <- ggplot(barley_site_year) +
    geom_point(aes(y = site, x = yield, color = year), size = 3)
p2 <- ggplot(barley_site_year) +
    geom_col(aes(x = site, y = yield, fill = year),
             size = 3,
             position = "dodge", width = .4) +
    coord_flip()
cowplot::plot_grid(p1, p2)

Because of the way we perceive bars, it is important to use a zero base line for bar charts.

Hair and Eye Color Data

A data set recording the distribution of hair and eye color and sex in 592 statistics students.

The data set is available as a cross-tabulation; as.data.frame converts it to a data frame.

HairEyeDF <- as.data.frame(HairEyeColor)
head(HairEyeDF)
##    Hair   Eye  Sex Freq
## 1 Black Brown Male   32
## 2 Brown Brown Male   53
## 3   Red Brown Male   10
## 4 Blond Brown Male    3
## 5 Black  Blue Male   11
## 6 Brown  Blue Male   50

Looking at the distribution of eye color:

eye <-
    group_by(HairEyeDF, Eye) %>%
    summarize(Freq = sum(Freq)) %>%
    ungroup()
ggplot(eye) +
    geom_col(aes(x = Eye, y = Freq), position = "dodge")

Mapping eye color to color in addition to the horizontal axis can help:

ggplot(eye) + geom_col(aes(x = Eye, y = Freq, fill = Eye), position = "dodge")

More sensible colors would be nice but requires a bit of work:

hazel_rgb <- col2rgb("brown") * 0.75 + col2rgb("green") * 0.25
hazel <- do.call(rgb, as.list(hazel_rgb / 255))

cols <- c(Blue = colorspace::lighten(colorspace::desaturate("blue", 0.3), 0.3),
          Green = colorspace::lighten("forestgreen", 0.1),
          Brown = colorspace::lighten("brown", 0.0001), ## 0.3?
          Hazel = colorspace::lighten(hazel, 0.3))

pb <- ggplot(eye) +
    geom_col(aes(x = Eye, y = Freq, fill = Eye), position = "dodge") +
    scale_fill_manual(values = cols)
pb

A stacked bar chart can also be useful:

psb <- ggplot(eye) +
    geom_col(aes(x = "", y = Freq, fill = Eye), color = "lightgrey") +
    scale_fill_manual(values = cols)
psb

A pie chart can be seen as a stacked bar chart in polar coordinates:

(pp <- psb + coord_polar("y"))

The axis and grid are not helpful; a theme adjustment can remove them:

(pp <- pp + theme_void())

Themes are a way to customize the non-data components of plots: i.e. titles, labels, fonts, background, gridlines, and legends. Themes can be used to give plots a consistent customized look.

The ggthemes package provides a number of themes to emulate the style of different publications, for example theme_wsj and theme_economist.

How well do bar charts and pie charts work?

cowplot::plot_grid(pb, pp)

Some questions:

• Which plot makes it easier to tell whether the proportion of brown-eyed students is larger or smaller that the proportion of blue-eyed students.

• Which plot makes it easier to tell whether these proportions are larger or smaller than 1/2 or 1/4 or 1/3?

Looking at the proportions within hair color and sex:

eye_hairsex <-
    group_by(HairEyeDF, Hair, Sex) %>%
    mutate(Prop = Freq / sum(Freq)) %>%
    ungroup()
p1 <- ggplot(eye_hairsex) +
    geom_col(aes(x = Eye, y = Prop, fill = Eye)) +
    scale_fill_manual(values = cols) +
    facet_grid(Hair~Sex)
p2 <- ggplot(eye_hairsex) +
    geom_col(aes(x = "", y = Prop, fill = Eye)) +
    scale_fill_manual(values = cols) +
    coord_polar("y")+facet_grid(Hair~Sex) +
    theme_void()
cowplot::plot_grid(p1, p2)

A more complete ggplot template:

ggplot(data = <DATA>) +
    <GEOM>(mapping = aes(<MAPPINGS>),
           stat = <STAT>,
           position = <POSITION>) +
    < ... MORE GEOMS ... > +
    <COORDINATE_ADJUSTMENT> +
    <SCALE_ADJUSTMENT> +
    <FACETING> +
    <THEME_ADJUSTMENT>

Perception and the Grammar of Graphics

river <- scan("data/river.dat")
rd <- data.frame(flow = river, month = seq_along(river))
(pp <- ggplot(rd) + geom_point(aes(x = month, y = flow)))

(pl <- ggplot(rd) + geom_line(aes(x = month, y = flow)))

pp + coord_fixed(3.5)

pl + coord_fixed(3.5)

A Simple Model of Visual Perception

The eyes acquire an image, which is processed through three stages of memory:

• Iconic memory
• Working memory, or short-term memory
• Long-term memory

The first processing stage of an image happens in iconic memory.

• Images remain in iconic memory for less than a second.
• Processing in iconic memory is massively parallel and automatic.
• This is called preattentive processing.

Preattentive processing is a fast recognition process.

Meaningful visual chunks are moved from iconic memory to short term memory.

• These chunks are used by conscious, or attentive, processing.
• Attentive processing often involves conscious comparisons or search.
• Short term memory is limited;
  • information is retained for only a few seconds;
  • only three or fours chunks can be held at a time.

Long term visual memory is built up over a lifetime, though infrequently used visual chunks may become lost.

Visual Design Implications

• Try to make as much use of preattentive features as possible.

• Recognize when preattentive features might mislead.

• For features that require attentive processing keep in mind that working memory is limited.

Some Terms for Describing Visualizations

• Data to be visualized contains variables or attributes measured on individual items or cases.

• Links are relationships that may exist among items, e.g. months within a year or countries within a continent.

• Marks are individual geometric entities used to represent items: points. bars, etc.

• Aesthetics or visual channels are the visual features of marks that can be used to encode attributes.

The aes(...) expressions establish the mapping between attributes and visual channels.

These ideas closely mirror the structure of the grammar of graphics as implemented in ggplot.

Munzner, T. (2014), Visualization Analysis and Design, CRC Press.

Wilkinson, L. (2005), The Grammar of Graphics, 2nd ed, Springer.

Channels and their Accuracy

A useful distinction among channels:

• Magnitude channels can reflect order and numeric values, e.g. position on an axis, length, area, brightness.

• Identity channels can distinguish different values but not reflect order, e.g. hue, shape, grouping.

Some channels are better at conveying information than others.

Munzner's ordering by accuracy:

Magnitude Channels (Ordered, Numerical)	Identity Channels (Categorical)
Position on common scale	Spatial grouping
Position on unaligned scale	Color hue
Length (1D size)	Shape
Tilt, angle
Area (2D size)
Depth (3D position)
Color luminance, saturation
Curvature, volume (3D size)

Line width is another channel; not sure there is agreement on its accuracy, but it is not high.

Visual Design Implications

Try to map the most important variables to the strongest channels.

Color

Color is very effective when used well.

But using color well is not easy.

Some of the issues:

• Perception depends on context.

• Simple color assignments may not separate equally well.

• Effectiveness may vary with the medium (screen, projector, print).

• Some people do not perceive the full specturm of colors.

• Grey scale printing.

• Some colors have cultural significance.

• Cultural significance may vary among cultures and with time.

Color perception is relative:

A note on rainbow colors.

Some tools for selecting palettes include:

• ColorBrewer; available in the RColorBrewer package.

• HCL Wizard; also available as hclwizard in the colorspace package.

A Grammar of Data Manipulation

The dplyr package provides a language, or grammar, for data manipulation.

The language contains a number of verbs that operate on tables.

The most commonly used verbs operate on a single data frame:

• select: pick variables by their names
• filter: choose rows that satisfy some criteria
• mutate: create transformed or derived variables
• arrange: reorder the rows
• summarize: collapse rows down to summaries

There are also a number of join verbs that merge several data frames into one.

The tidyr package provides additional verbs, such as pivot_longer and pivot_wider for reshaping data frames.

The single table verbs can also be used with group_by to work a group at a time instead of applying to the entire data frame.

The design of dplyr is strongly motivated by SQL.

More Examples

These examples start with raw data as you might receive it from a researcher, and involve reading and cleaning the data.

Wind Turbines in Iowa

There are many wind turbines in Iowa. Data is available from the U.S. Wind Turbine Database. A snapshot is available is here as a CSV file.

CSV files are a common form of data exchange.

• They are simple text files that are intended to be written and read by a computer.

• Some CSV files include a header and a footer that need to he handled.

• One issue is that a comma isn't a good separator in countries where it is the decimal separator!

• A CSV file can be read using read.csv or readr::read_csv.

Reading the wind turbine data:

wind_turbines <- read.csv("data/us_wind.csv", comment = "#")

Some data cleaning is needed.

Focus on the wind turbines in IOWA (19 is the FIPS county code for Iowa):

wt_IA <- filter(wind_turbines, t_fips %/% 1000 == 19)

Drop entries with missing longitude or latitude values:

wt_IA <- filter(wt_IA, ! is.na(xlong), ! is.na(ylat))

Some missing year values are encoded as -9999; replace these with NA:

wt_IA <- mutate(wt_IA, p_year = replace(p_year, p_year < 0, NA))

To show the locations of wind turbines on a map, load some map data:

iowa_sf <- sf::st_as_sf(maps::map("county", "iowa", plot = FALSE, fill = TRUE))

To show the locations of wind turbines on a map, load some map data:

iowa_sf <-
    sf::st_as_sf(maps::map("county", "iowa",
                           plot = FALSE,
                           fill = TRUE))
p <- ggplot() +
    geom_sf(data = iowa_sf) +
    ggthemes::theme_map()
p

Locations for all wind turbines in iowa:

p + geom_point(aes(xlong, ylat), data = wt_IA)

Using color to show when the wind turbines were built:

year_brk <-c(0, 2005, 2010, 2015, 2020)
year_lab <- c("before 2005",
              "2005-2009",
              "2010-2014",
              "2015-2020")
wt_IA <-
    mutate(wt_IA,
           year = cut(p_year,
                      breaks = year_brk,
                      labels = year_lab,
                      right = FALSE))
p + geom_point(aes(xlong,
                   ylat,
                   color = year),
               data = wt_IA,
               size = 3)

Cancer Map

The website http://www.cancer-rates.info/ia provides data on cancer incidence for a number of different cancers in Iowa. The data for lung and bronchus cancer in 2011 are available in a csv file in the project.

We can read the file with read_csv from the readr package.

Looking at the file shows some things that need to be cleaned up:

• Two header lines at the beginning
• Some footer lines.
• Some values codes as ~.

The header can be handled by using skip = 2 in the read_csv call:

fname <- "data/Invasive-Cancer-Incidence-Rates-by-County-in-Iowa-Lung-and-Bronchus-2011.csv"
d <- read_csv(fname, skip = 2)
## Warning: 3 parsing failures.
## row col  expected    actual                                                                                file
## 101  -- 7 columns 1 columns 'data/Invasive-Cancer-Incidence-Rates-by-County-in-Iowa-Lung-and-Bronchus-2011.csv'
## 102  -- 7 columns 1 columns 'data/Invasive-Cancer-Incidence-Rates-by-County-in-Iowa-Lung-and-Bronchus-2011.csv'
## 103  -- 7 columns 1 columns 'data/Invasive-Cancer-Incidence-Rates-by-County-in-Iowa-Lung-and-Bronchus-2011.csv'
head(d)
## # A tibble: 6 x 7
##   County  `Population at … Cases `Crude Rate` `Age-adjusted R… `95% Confidence …
##   <chr>              <dbl> <chr> <chr>        <chr>            <chr>            
## 1 Union              12570 20    159.11       115.82           69.86            
## 2 Ringgo…             5098 11    215.77       115.03           54.48            
## 3 Monroe              8044 12    149.18       109.99           56.24            
## 4 Page               15926 25    156.98       109.20           70.12            
## 5 Montgo…            10655 17    159.55       99.45            57.41            
## 6 Adams               3996 6     150.15       95.84            35.14            
## # … with 1 more variable: 95% Confidence Interval-Upper Limit <chr>

Let's focus on a few variables and give them more convenient names:

d <- select(d, county = 1, population = 2, count = 3, crude_rate = 4)

The footer needs to be removed:

tail(d)
## # A tibble: 6 x 4
##   county                                             population count crude_rate
##   <chr>                                                   <dbl> <chr> <chr>     
## 1 Butler                                                  14960 5     33.42     
## 2 Winneshiek                                              21045 ~     ~         
## 3 STATE                                                 3065223 2368  77.25     
## 4 Note: All rates are per 100,000. Rates are age-ad…         NA <NA>  <NA>      
## 5 Rates generated on Jun 12, 2019.                           NA <NA>  <NA>      
## 6 Based on data released Nov 2017.                           NA <NA>  <NA>

One way to remove the footer:

d <- filter(d, ! is.na(population))
d <- filter(d, county != "STATE")
tail(d)
## # A tibble: 6 x 4
##   county     population count crude_rate
##   <chr>           <dbl> <chr> <chr>     
## 1 Lyon            11720 ~     ~         
## 2 Kossuth         15392 7     45.48     
## 3 Palo Alto        9360 ~     ~         
## 4 Grundy          12474 5     40.08     
## 5 Butler          14960 5     33.42     
## 6 Winneshiek      21045 ~     ~

Changing count and crude_rate to numeric changes the ~ entries to missing values (NA) values:

d <- mutate(d, count = as.numeric(count), crude_rate = as.numeric(crude_rate))
## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion

## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion

In this case there are no zero case values; two ways to check:

count(d, count == 0)
## # A tibble: 2 x 2
##   `count == 0`     n
##   <lgl>        <int>
## 1 FALSE           95
## 2 NA               4
any(d$count == 0, na.rm = TRUE)
## [1] FALSE

It might be reasonable to assume these values where zero, so replace them with zeros:

d <- replace_na(d, list(count = 0, crude_rate = 0))

A choropleth map uses color or shading to represent values measured for different geographic regions.

We will need to merge, or left join, the cancer data with the map date we loaded for the wind turbine map.

For Iowa this can be done with the county name, but some care is needed.

d$county[1]
## [1] "Union"
iowa_sf$ID[1]
## [1] "iowa,adair"

d <- mutate(d, cname = county, county = tolower(county))
iowa_sf <- mutate(iowa_sf, county = sub("iowa,", "", ID))

setdiff(d$county, iowa_sf$county)
## [1] "o'brien"
setdiff(iowa_sf$county, d$county)
## [1] "obrien"

d <- mutate(d, county = sub("'", "", county))

setdiff(d$county, iowa_sf$county)
## character(0)
setdiff(iowa_sf$county, d$county)
## character(0)

Define rate1K variable as the number of cases per 1000 inhabitants and left join the data to the polygons:

d <- mutate(d, rate1K = 1000 * (count / population))
md <- left_join(iowa_sf, d, "county")
head(md)
## Simple feature collection with 6 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -95.10526 ymin: 40.60552 xmax: -91.06018 ymax: 43.51041
## CRS:           EPSG:4326
##               ID    county population count crude_rate     cname    rate1K
## 1     iowa,adair     adair       7565    10     132.19     Adair 1.3218771
## 2     iowa,adams     adams       3996     6     150.15     Adams 1.5015015
## 3 iowa,allamakee allamakee      14204     7      49.28 Allamakee 0.4928189
## 4 iowa,appanoose appanoose      12863    11      85.52 Appanoose 0.8551660
## 5   iowa,audubon   audubon       6019     5      83.07   Audubon 0.8307028
## 6    iowa,benton    benton      26121    28     107.19    Benton 1.0719345
##                             geom
## 1 MULTIPOLYGON (((-94.24583 4...
## 2 MULTIPOLYGON (((-94.70992 4...
## 3 MULTIPOLYGON (((-91.22634 4...
## 4 MULTIPOLYGON (((-92.63009 4...
## 5 MULTIPOLYGON (((-95.10526 4...
## 6 MULTIPOLYGON (((-92.06286 4...

A simple map:

library(ggthemes)
library(viridis)
## Loading required package: viridisLite
ggplot(md) + geom_sf(aes(fill = rate1K))

An improved version:

library(ggthemes)
library(viridis)
ggplot(md) +
    geom_sf(aes(fill = rate1K),
            color = "grey") +
    scale_fill_viridis(name = "Rate per 1000") +
    theme_map()

A simple interactive version using plotly:

mdl <- mutate(md,
              label = paste(cname, round(rate1K, 1), population, sep = "\n"))
p <- ggplot(mdl) +
    geom_sf(aes(fill = rate1K,
                text = label), 
            color = "grey") +
    scale_fill_viridis(name = "Rate per 1000") +
    theme_map()
## Warning: Ignoring unknown aesthetics: text

plotly::ggplotly(p, tooltip = "text")

The leaflet package supports more sophisticated interactive maps:

library(leaflet)
pal <- colorNumeric(
    palette = "viridis",
    domain = md$rate1K)
lab <- lapply(paste0(tools::toTitleCase(md$county), "<BR>",
                     "Rate: ", round(md$rate1K, 1), "<BR>",
                     "Pop: ", scales::comma(md$population,
                                            accuracy = 1)),
              htmltools::HTML)
leaflet(sf::st_transform(md, 4326)) %>%
    addPolygons(weight = 2,
                color = "grey",
                fillColor = ~ pal(rate1K),
                fillOpacity = 1,
                highlightOptions =
                    highlightOptions(color = "white",
                                     weight = 2,
                                     bringToFront = TRUE),
                label = lab) %>%
    addLegend(pal = pal, values = ~ rate1K)

Unemployment Map

Local Area Unemployment Statistics page from the Bureau of Labor Statistics makes available county-level monthly unemployment data for a 14-month window. The file for February 2020 through March 2021 is available is available at http://www.stat.uiowa.edu/~luke/data/laus/laucntycur14-2020.txt and in the project data folder.

This file is a text file but uses a non-standard separator. It is designed for human readability and uses a comma as a thousands separator or grouping mark. It also includes header and footer information. It is still reasonably easy to read in.

One way to read the data into R is:

lausURL <- "data/laucntycur14-2020.txt"
lausUS <- read.table(lausURL,
                     col.names = c("LAUSAreaCode", "State", "County",
                                   "Title", "Period",
                                   "LaborForce", "Employed",
                                   "Unemployed", "UnempRate"),
                     quote = '"', sep = "|", skip = 6,
                     stringsAsFactors = FALSE, strip.white = TRUE,
                     fill = TRUE)
footstart <- grep("------", lausUS$LAUSAreaCode)
lausUS <- lausUS[1:(footstart - 1),]

It may be useful to be able to access the county name and state name separately:

lausUS <- separate(lausUS, Title, c("cname", "scode"),
                   sep = ", ", fill = "right")

The UnempRate variable is read as character data because of missing value encoding, so needs to be converted to numeric:

lausUS <- mutate(lausUS, UnempRate = as.numeric(UnempRate))
## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion

Check for missing values:

select_if(lausUS, anyNA) %>% names()
## [1] "scode"     "UnempRate"

The state code is missing for the District of Columbia:

select(lausUS, cname, scode) %>%
    filter(is.na(scode)) %>%
    unique()
##                  cname scode
## 1 District of Columbia  <NA>

March and April 2020 numbers were not available for Puerto Rico:

select(lausUS, scode, Period, UnempRate) %>%
    filter(is.na(UnempRate)) %>%
    unique()
##    scode Period UnempRate
## 1     PR Mar-20        NA
## 79    PR Apr-20        NA

To compute the national monthly unemployment rates over this period we need some more data cleaning:

lausUS <- mutate(lausUS,
                 Period = fct_inorder(Period),
                 LaborForce = as.numeric(gsub(",", "", LaborForce)),
                 Unemployed = as.numeric(gsub(",", "", Unemployed)))
## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion

## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion

Unemployment during this period was affected significantly by the COVID-19 pandemic. A plot shows a large spike in April 2020:

group_by(lausUS, Period) %>%
    summarize(Unemployed = sum(Unemployed, na.rm = TRUE),
              LaborForce = sum(LaborForce, na.rm = TRUE),
              UnempRate = 100 * (Unemployed / LaborForce)) %>%
    ggplot(aes(Period, UnempRate, group = 1)) +
    geom_line()

A choropleth map can be used to look at how the impact was distributed across the country.

To show unemployment rates on a map we need to merge the unemployment data with map data.

To match county unemployment data and county shape data it is safer to use the numeric FIPS county code. This can be added with

lausUS <- mutate(lausUS, fips = State * 1000 + County)

Shape data for US counties can be obtained from a number of sources in a number of different formats.

Here is one approach:

counties_sf <- sf::st_as_sf(maps::map("county", plot = FALSE, fill = TRUE))
county.fips <-
    mutate(maps::county.fips, polyname = sub(":.*", "", polyname)) %>%
    unique()
counties_sf <- left_join(counties_sf, county.fips, c("ID" = "polyname"))
states_sf <-  sf::st_as_sf(maps::map("state", plot = FALSE, fill = TRUE))

Some summaries over the period can be computed as

summaryUS <- group_by(lausUS, County, State, fips) %>%
    summarize(avg_unemp = mean(UnempRate, na.rm = TRUE),
              max_unemp = max(UnempRate, na.rm = TRUE),
              apr_unemp = UnempRate[Period == "Apr-20"]) %>%
    ungroup()
## `summarise()` has grouped output by 'County', 'State'. You can override using the `.groups` argument.
head(summaryUS)
## # A tibble: 6 x 6
##   County State  fips avg_unemp max_unemp apr_unemp
##    <int> <int> <dbl>     <dbl>     <dbl>     <dbl>
## 1      1     1  1001      4.65      10.9      10.9
## 2      1     4  4001     12.7       17.6      16.6
## 3      1     5  5001      4.06       5.2       4.8
## 4      1     6  6001      8.8       14.6      14.6
## 5      1     8  8001      8.3       12.7      12.7
## 6      1     9  9001      8.31      11.7       8.2

A choropleth map of the April 2020 unemployment rates:

left_join(counties_sf, summaryUS, "fips") %>%
    ggplot() +
    geom_sf(aes(fill = apr_unemp)) +
    scale_fill_viridis(name = "Rate", na.value = "red") +
    theme_map() +
    geom_sf(data = states_sf, col = "grey", fill = NA)

Using a very visible color for missing data is useful, at least during exploration.

anti_join can show the county geometry that does not have an entry in the unemployment data:

anti_join(counties_sf, summaryUS, "fips")
## Simple feature collection with 1 feature and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -103.0121 ymin: 42.99475 xmax: -102.0782 ymax: 43.68803
## CRS:           EPSG:4326
##                     ID  fips                           geom
## 1 south dakota,shannon 46113 MULTIPOLYGON (((-102.8115 4...

Shannon County, SD (FIPS 46113), was renamed to Oglala Lakota County in June 2015 and given a new FIPS code, 46102.

The geometry data table needs to be updated:

counties_sf <- mutate(counties_sf, fips = replace(fips, fips == 46113, 46102))

With the updated data the map is now complete:

left_join(counties_sf, summaryUS, "fips") %>%
    ggplot() +
    geom_sf(aes(fill = apr_unemp)) +
    scale_fill_viridis(name = "Rate", na.value = "red") +
    theme_map() +
    geom_sf(data = states_sf, col = "grey", fill = NA)

Gapminder Childhood Mortality Data

The gapminder package provides a subset of the data from the Gapminder web site. Additional data sets are available.

• A data set on childhood mortality is available locally as a csv file or an Excel file. The Excel file is also available in the project data folder.

• The numbers represent number of deaths within the first five years per 1000 births.

Many researchers like to manage their data in a spreadsheet. Being able to read such a sheet directly greatly helps keeping the workflow reproducible.

Many spreadsheets contain header, footers, and other annotations to aid a human viewer.

As long as the data are in a rectangular region it is usually not hard to extract them programmatically.

Loading the data:

library(readxl)
gcm <- read_excel("data/gapminder-under5mortality.xlsx")
names(gcm)[1]
## [1] "Under five mortality"
names(gcm)[1] <- "country"
head(gcm)
## # A tibble: 6 x 217
##   country         `1800.0` `1801.0` `1802.0` `1803.0` `1804.0` `1805.0` `1806.0`
##   <chr>              <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
## 1 Abkhazia             NA       NA       NA       NA       NA       NA       NA 
## 2 Afghanistan         469.     469.     469.     469.     469.     469.     470.
## 3 Akrotiri and D…      NA       NA       NA       NA       NA       NA       NA 
## 4 Albania             375.     375.     375.     375.     375.     375.     375.
## 5 Algeria             460.     460.     460.     460.     460.     460.     460.
## 6 American Samoa       NA       NA       NA       NA       NA       NA       NA 
## # … with 209 more variables: 1807.0 <dbl>, 1808.0 <dbl>, 1809.0 <dbl>,
## #   1810.0 <dbl>, 1811.0 <dbl>, 1812.0 <dbl>, 1813.0 <dbl>, 1814.0 <dbl>,
## #   1815.0 <dbl>, 1816.0 <dbl>, 1817.0 <dbl>, 1818.0 <dbl>, 1819.0 <dbl>,
## #   1820.0 <dbl>, 1821.0 <dbl>, 1822.0 <dbl>, 1823.0 <dbl>, 1824.0 <dbl>,
## #   1825.0 <dbl>, 1826.0 <dbl>, 1827.0 <dbl>, 1828.0 <dbl>, 1829.0 <dbl>,
## #   1830.0 <dbl>, 1831.0 <dbl>, 1832.0 <dbl>, 1833.0 <dbl>, 1834.0 <dbl>,
## #   1835.0 <dbl>, 1836.0 <dbl>, 1837.0 <dbl>, 1838.0 <dbl>, 1839.0 <dbl>,
## #   1840.0 <dbl>, 1841.0 <dbl>, 1842.0 <dbl>, 1843.0 <dbl>, 1844.0 <dbl>,
## #   1845.0 <dbl>, 1846.0 <dbl>, 1847.0 <dbl>, 1848.0 <dbl>, 1849.0 <dbl>,
## #   1850.0 <dbl>, 1851.0 <dbl>, 1852.0 <dbl>, 1853.0 <dbl>, 1854.0 <dbl>,
## #   1855.0 <dbl>, 1856.0 <dbl>, 1857.0 <dbl>, 1858.0 <dbl>, 1859.0 <dbl>,
## #   1860.0 <dbl>, 1861.0 <dbl>, 1862.0 <dbl>, 1863.0 <dbl>, 1864.0 <dbl>,
## #   1865.0 <dbl>, 1866.0 <dbl>, 1867.0 <dbl>, 1868.0 <dbl>, 1869.0 <dbl>,
## #   1870.0 <dbl>, 1871.0 <dbl>, 1872.0 <dbl>, 1873.0 <dbl>, 1874.0 <dbl>,
## #   1875.0 <dbl>, 1876.0 <dbl>, 1877.0 <dbl>, 1878.0 <dbl>, 1879.0 <dbl>,
## #   1880.0 <dbl>, 1881.0 <dbl>, 1882.0 <dbl>, 1883.0 <dbl>, 1884.0 <dbl>,
## #   1885.0 <dbl>, 1886.0 <dbl>, 1887.0 <dbl>, 1888.0 <dbl>, 1889.0 <dbl>,
## #   1890.0 <dbl>, 1891.0 <dbl>, 1892.0 <dbl>, 1893.0 <dbl>, 1894.0 <dbl>,
## #   1895.0 <dbl>, 1896.0 <dbl>, 1897.0 <dbl>, 1898.0 <dbl>, 1899.0 <dbl>,
## #   1900.0 <dbl>, 1901.0 <dbl>, 1902.0 <dbl>, 1903.0 <dbl>, 1904.0 <dbl>,
## #   1905.0 <dbl>, 1906.0 <dbl>, …

This data set is in wide format.

A long version is useful for working with ggplot.

tgcm <- pivot_longer(gcm, -1, names_to = "year", values_to = "u5mort")
head(tgcm)
## # A tibble: 6 x 3
##   country  year   u5mort
##   <chr>    <chr>   <dbl>
## 1 Abkhazia 1800.0     NA
## 2 Abkhazia 1801.0     NA
## 3 Abkhazia 1802.0     NA
## 4 Abkhazia 1803.0     NA
## 5 Abkhazia 1804.0     NA
## 6 Abkhazia 1805.0     NA
tgcm <- mutate(tgcm, year = as.numeric(year))
head(tgcm)
## # A tibble: 6 x 3
##   country   year u5mort
##   <chr>    <dbl>  <dbl>
## 1 Abkhazia  1800     NA
## 2 Abkhazia  1801     NA
## 3 Abkhazia  1802     NA
## 4 Abkhazia  1803     NA
## 5 Abkhazia  1804     NA
## 6 Abkhazia  1805     NA

Some explorations:

p <- ggplot(tgcm) +
    geom_line(aes(year, u5mort, group = country), alpha = 0.3)
p
## Warning: Removed 18644 row(s) containing missing values (geom_path).

plotly::ggplotly(p)

Some selected countries:

countries <- c("United States", "United Kingdom", "Germany", "China", "Egypt")
filter(tgcm, country %in% countries) %>%
    ggplot() +
    geom_line(aes(x = year, y = u5mort, color = country))

Examining the missing values:

tgcm_miss <-
    group_by(tgcm, country) %>%
    summarize(anyNA = any(is.na(u5mort))) %>%
    filter(anyNA) %>%
    pull(country)

p <- ggplot(filter(tgcm, country %in% tgcm_miss)) +
    geom_line(aes(x = year, y = u5mort, group = country), na.rm = TRUE)
p

plotly::ggplotly(p)

This work is licensed under the CC BY-NC 3.0 Creative Commons License.